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ANNOUNCEMENT 


The  Plasma  Theory  and  Simulation  6roup  in  the  EECS  Department 
and  ERL  is  pleased  to  announce  receipt  of  three  awards  to  its  members. 

The  first  is  the  appointment  of  Dr.  Liu  Chen  as  Visiting  Miller  Research 
Professor  for  an  eight  week  visit  during  the  1987-88  academic  year,  to 
be  hosted  by  our  group  in  EECS,  with  Prof.  Birdsall  as  faculty  sponsor.  Dr. 
Chen  comes  to  us  from  Princeton  University  where  he  is  Deputy  Head  of  the 
theory  group  at  the  Princeton  Plasma  Physics  Laboratory.  Dr  Chen  received 
his  Ph.D.  while  in  our  group  in  the  early  1970’s.  It  is  expected  that  he  will 
arrive  in  the  Fall  1987  and  interact  with  the  plasma  community  which 
reaches  across  several  department  lines,  as  well  as  with  researchers  at 
LBL  and  LLNL. 

The  second  is  the  award  to  Lou  Ann  Schwager  of  a  U.S.  Dept,  of  Energy 
Fusion  Energy  Postdoctoral  Award,  with  appointment  to  be  at  the  Sandia 
National  Laboratories,  Livermore.  Ms.  Schwager  is  now  completing  her 
dissertation  on  plasma-sheath-wall  dynamics,  Including  secondary 
electron  emission  and  reflected  ions,  working  in  our  PTSG,  a  student  in 
NE.  She  will  take  up  her  post-doc  about  September  1.  It  is  expected  that 
she  will  split  her  time  between  Sandia,  where  Drs.  Walter  Bauer  and  Ken 
Wilson  have  excellent  plasma-surface  experimental  facilities,  and  PTSG, 
where  she  will  continue  with  her  simulations.  The  award  is  one  of  three 
made  nationally. 

The  third  is  the  award  to  Richard  Procassini  of  a  DAAD  (Deutscher 
Akademischer  Austauschdienst)  scholarship  to  continue  his  graduate 
plasma  studies  for  a  year  at  the  University  of  Munich,  Max  Planck  Involute 
for  Plasma  Physics,  at  Garching,  working  with  our  friend  Dr.  Roland 
Chodura.  Mr.  Procassini  is  an  NE  graduate  student  in  PTSG,  with  much  of 
his  research  supported  by  an  IR&D  grant  from  LLNL.  He  is  working  closely 
with  Dr.  Bruce  Cohen  at  LLNL,  on  a  large  (full  fusion  machine  particle 
simulation)  code  TES5,  developed  by  Dr.  Cohen  and  Prof.  Birdsall.  There  is 
similar  work  at  MPIPP,  especially  on  the  addition  of  Coulomb  collisions. 

The  Plasma  Theory  and  Simulation  Group  is  now  completing  its  20th 
year,  dating  back  to  Prof.  Birdsall's  return  from  sabbatical  in  Japan  with 
Prof.  A.  Hasegawa  and  T.  Kamimura,  and  the  appointment  of  Dr.  A.  B. 

Langdon  to  be  our  first  post-doc.  Fourteen  Ph.D.s  and  six  Masters  have 
been  granted  in  that  period  to  students  in  the  Group,  and  there  have  been 
six  post-docs  with  two  and  three  year  terms.  The  major  text  in  the  area, 
Plasma  Physics  via  Computer  Simulation  (1985),  by  Birdsall  and  Langdon, 
has  included  much  of  the  work  of  the  Group,  in  methods  and  in  applications. 
Students  and  post-docs  have  come  from  EECS,  NE,  Physics,  Math  and 
Astronomy  and  have  moved  on  to  work  in  a  wide  variety  of  engineering  and 
scientific  activities. 

Looking  forward  to  a  second  20  years! 

20  April  1987 


Charles  K.  (Ned)  Birdsall 

With  apologies  for  the  gild  on  the  Illy! 
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C .  Implicit  particle  simulation  of  velocity  space  transport  and 
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SECTION  b  GENERAL  PLASMA  THEORY  AND  SIMULATION 


A.  Computer  Simulation  of  Bounded  Plasma  Systems 

William  S.  Lawson 

Abstract 

The  physical  and  numerical  problems  of  kinetic  simulation  of  a  bounded  electrostatic 
plasma  system  in  one  planar  dimension  are  examined,  and  solutions  to  them  are  presented. 
These  problems  include  particle  absorption,  reflection  and  emission  at  boundaries,  the  solu¬ 
tion  of  Poisson's  equation  under  non-periodic  boundary  conditions,  and  the  treatment  of 
an  external  circuit  connecting  the  boundaries.  Some  comments  are  also  made  regarding  the 
problems  of  higher  dimensions.  The  methods  which  are  described  here  are  implemented  in 
a  code  named  PDW1.  which  is  available  from  Professor  C.  K.  Birdsall.  Plasma  Theory  and 
Simulation  Group,  Cory  Hall.  University  of  California.  Berkeley.  CA  94720. 

This  work  is  to  be  issued  as  Memorandum  No.  UCB/ERL  M87/14.  March  5. 1987. 

B.  Artificial  Cooling  of  Trapped  Electrons  in  Bounded  Simulations 

William  S.  Lawson 

Abstract 

An  explanation  is  proposed  for  an  artificial  cooling  effect  seen  in  electrostatic 
particle-in-cell  plasma  simulations.  Further  simulations  are  done  which  test  and  support 
the  explanation. 

A  report  with  the  title  shown  above  is  in  preparation.  Memorandum  No.  UCB/ERL 
M87/34.  as  well  as  for  journal  submission.  The  final  title  may  be  altered  slightly. 


Projects  reported  as  A.B.C.D.E.F  in  the  previous  QPR  (HI  &  IV.  1985)  have  been  com¬ 
pleted  and  final  reports  have  been  issued. 
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SECTION  m  PLASMA-WALL  PHYSICS,  THEORY  AND  SIMULATION 


A.  Collector  Sheath,  Presheath,  and  Source  Sheath  in  a  Collisionless,  Finite  Ion 

Temperature  Plasma 

Lou  Ann  Schwager 

B.  Potential  Drop  and  Transport  in  a  Bounded  Plasma  with  Ion  Reflection 

at  the  Collector 

Lou  Ann  Schwager 

C.  Potential  Drop  and  Transport  in  a  Bounded  Plasma  with  Secondary  Electron 

Emission  at  the  Collector 

Lou  Ann  Schwager 

Final  reports  with  the  above  titles  are  in  preparation  for  ERL  memos  and  journal 
submissions.  The  final  titles  may  be  altered  slightly. 

D.  Magnetized  Sheath;  Kelvin-Helmholtz  Instability;  Long-lived  Vortices 

Dr.  Kim  Theilhaber 

In  performing  two-dimensional  simulations  of  a  magnetized  sheath,  we  expect  to 
reveal  important  physical  behavior  which  is  not  described  by  the  one-dimensional  simula¬ 
tions.  in  particular  time-varying  electric  field  fluctuations  sustained  by  plasma  instabili¬ 
ties.  and  intrinsically  two-dimensional  particle  transport  across  the  discharge.  We  have 
been  using  our  two-dimensional  particle  simulation  code  ES2  to  model  a  bounded  plasma 
in  a  plane  perpendicular  to  the  magnetic  field,  along  the  midsection  shown  in  Fig.  1.  page  8. 
QPR  III  &  IV.  1985.  In  this  configuration,  the  walls  of  the  device  lie  parallel  to  the  y  direc¬ 
tion.  and  are  separated  by  a  distance  in  the  x  direction,  with  the  magnetic  field  parallel  to 
z.  Because  the  electric  fields  in  the  plasma  are  spatially  non-uniform,  and  give  rise  to  a 
highly  sheared  E  x  B  drift  of  electrons  and  ions  in  the  y-direction.  the  discharge  will  be 
unstable  to  the  so-called  "Kelvin-Helmholtz”  instability,  an  instability  which  feeds  on  the 
free  energy  contained  in  the  sheared  flow.  In  the  long  term,  starting  from  an  initially  lam¬ 
inar  flow,  the  Kelvin-Helmholtz  instability  drives  the  particle  motion  into  a  set  of  vor¬ 
tices.  through  which  some  but  not  all  of  the  free  energy  in  the  shear  layer  is  dissipated. 

A  snapshot  from  an  ES2  simulation,  in  which  the  Kelvin-Helmholtz  instability  was 
present,  is  shown  in  Fig.  1.  This  figure  shows  the  electrostatic  potential  contour  lines  for 
the  final,  quasi-steady  state  of  the  instability,  at  a>cit  =400.  in  a  system  which  was  ini¬ 
tially  filled  with  a  uniform  plasma..  This  particular  situation  did  not  faithfully  reproduce 
the  conditions  of  the  magnetized  sheath,  in  particular  because  there  was  no  external  bias 
imposed  across  the  system;  rather,  it  was  meant  as  a  preliminary  exploration  of  the  model. 
In  the  simulation,  the  wall  potentials  were  left  floating,  and  the  initial  shear  flow  was  pro¬ 
vided  by  the  electric  field  generated  by  charge  deposition  at  the  walls,  due  to  the  depletion, 
over  a  layer  of  the  order  of  an  ion  gyro  radius  in  width,  of  ions  near  the  wall.  The  system 
pictured  in  Fig.  1  is  about  5  ion  gyro  radii  across,  and  about  20  ion  gyro  radii  long  in  the 
direction  parallel  to  the  walls.  While  the  initial  vortices,  which  develop  out  of  the 
Kelvin-Helmholtz  instability  at  otcit  =*20  are  localized  near  the  walls,  in  Fig.  1  we  are 
looking  at  the  potential  at  the  later  time  utc,t  =400.  and  the  vortices  have  expanded  in  x  so 
as  to  occupy  the  entire  system.  The  potential  perturbations  induced  by  the  vortices  are 
sizeable.  e<fr/  T,  =0.5.  They  carry  along  comparable  perturbations  in  the  particle  density, 
within  /  n0  =0.5. 


The  most  striking  feature  of  the  vortices  is  their  great  stability:  the  vortices  on  the 
left  move  with  a  steady  motion  parallel  to  the  wall,  in  accordance  with  the  EiB  drift, 
while  those  near  the  right  wall  stream  in  the  opposite  direction  with  equal  speed.  The 
vortices  overlap  and  pass  through  each  other  with  little  effect  on  their  shapes  and  ampli¬ 
tudes. 

An  important  consequence  of  the  presence  of  vortex  motion  in  the  plasma  is  the  rate 
of  particle  transport  to  the  walls,  which  is  enhanced  over  the  collisional  rate.  In  the  simu¬ 
lation  of  Fig.  1.  the  enhanced  transport  resulted  in  the  depletion  of  almost  all  particles  in 
the  simulation  region,  over  a  time  scale  much  shorter  than  the  collisional  time  scale.  We 
are  presently  studying  the  dependence  of  this  "anomalous"  transport  on  the  vortex  struc¬ 
ture  and  on  the  plasma  parameters.  Another  important  feature  of  the  vortices,  which 
should  be  observable  in  an  experiment,  is  that  they  produce  a  high-frequency  power  spec¬ 
trum.  which  peaks  at  v0 .  where  vD  is  the  characteristic  drift  velocity,  and  ky  the 
characteristic  wave  vector  of  the  potential  fluctuations  parallel  to  the  wall. 

We  expect  behavior  similar  to  that  discussed  above  in  a  fuller  model  of  the  magnet¬ 
ized  sheath,  with  nonzero  bias.  The  enhanced  transport  due  to  the  vortices  should  be  espe¬ 
cially  important  in  the  magnetized  sheath,  as  we  believe  that  this  might  modify  the 
sputtering  efficiency  of  the  device.  Two-dimensional  particle  simulations  are  now  under 
way  to  study  how  this  efficiency  might  be  affected,  and  how  it  depends  on  the  plasma 
parameters  such  as  density  and  magnetic  field. 


fit  ■  ( X)  :  Contour  fist  of  tbo  electroatatlc  potential 
ep(s.y),  at  W„t»  100,  In  a  bounded  plaene  (Inula don, 
oelng  the  code  "XSl".  The  confining  walla  are  at  s»0  and 
said,  the  asternal  aagnetlc  field  la  perpendicular  to  the 
plane  of  the  emulation  region.  The  left-hand  aortas  la 
etraaalng  down  the  wall  with  a  velocity  of  about  0.1  v  . 


SECTION  m:  CODE  DEVELOPMENT  AND  SOFTWARE  DISTRIBUTION 


A.  ES2  User's  Manual— Version  1 

Dr.  Kim  Theilhaber 

Abstract 

The  two-dimensional,  electrostatic  plasma  simulation  code  “ES2” 
is  described.  This  includes  a  synopsis  of  the  physical  basis  of  the 
model  and  a- summary  of  the  numerical  schemes  used  to  implement 
it.  This  is  followed  by  two  examples  of  applications,  one  showing  the 
.  formation  of  unmagnetized  double  layers,  the  other  illustrating  the 
formation  of  vortices  in  a  configuration  perpendicular  to  the  magnetic 
field.  Finally,  the  main  section  of  the  manual  is  presented:  a  copy  of 
the  “on-line”  manual,  in  which  a  detailed  description  of  the  various 
options  and  parameters  entering  into  ES2  is  given. 

This  work  has  been  issued  as  an  ERL  Memo  No.  UCB/ERL  M87/23,  May  11. 1987. 

B.  Addition  of  Electron-Neutral  Elastic  Collisions 

Perry  Gray  (write-up  by  Prof.  C.  K.  Birdsall) 

Electron-neutral  collisions  may  be  accounted  for  in  simulations,  as  for  application  to 
low-pressure  discharges.  Let  the  electrons  collide  with  hard  sphere  neutral  atoms,  so  that 
the  collisions  are  elastic  (no  momentum  transfer)  and  the  electrons  are  scattered  randomly 
with  a  uniform  distribution  in  velocity  angle,  independent  of  their  initial  velocity  angle. 

The  details  of  this  program  and  the  checks  performed  will  be  given  in  a  later  QPR. 


C .  Implicit  Particle  Simulation  of  Velocity  Space 
Transport  and  Self-Consistent  Electric 
Potentials  in  Magnetized  Plasmas 

TESS  Code 

R.  J.  Procasaini,  C.  K.  Birdaoll  and  E.  C.  Morse 

and 

B.  I.  Cohen  (Lawrence  Livermore  National  Laboratory) 
Introduction 

We  are  developing  a  one-spatial-dimension,  three-velocity-component  implicit 
particle  simulation  code  incorporating  binary  particle  Coulomb  collisions  self-con 
sistently  to  simulate  velocity  space  transport  and  electric  potentials  in  magnetized 
plasmas.  The  simulation  model  addresses  long  transport  time  scales,  which  are  not 
accessible  to  explicit  algorithms,  and  retains  essential  kinetic  details.  This  new  sim¬ 
ulation  tool  is  an  efficient  alternative  to  Fokker-Planck  codes.  The  implicit  potential 
solver,  particle  pusher,  and  collision  routines  have  been  successfully  implemented 
and  tested.  Rudimentary  simulations  of  a  Q-machine  and  a  simple  tandem  mirror 
configuration  have  been  performed.  In  the  near  future,  we  will  add  physics  packages 
to  model  radiofrequency  wave  heating,  plasma  fueling  by  gas-puffing  ,  neutral  beam 
heating  and  fueling,  and  particle  losses.  At  that  point,  we  will  be  able  to  simulate 
thermal  barrier  tandem  mirror  operation.  Figure  1  shows  a  schematic  representa¬ 
tion  of  the  timestep  cycle  of  the  code,  including  the  particle  pusher,  field  solver, 
and  collision,  heating  and  particle  loss  routines. 

Objectives 

We  are  combining  recently  introduced  implicit  particle  simulation  methods1 
with  Monte  Carlo  model  of  binary  particle  collisions,  neutral  beam  ionization  and 
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Figure  1 .  The  TESS  Computational  Sequence. 


-  6  - 


charge  exchange,  and  radiofrequency  wave  heating  to  perform  simulations  of  both 
velocity-space  transport  and  the  evolution  of  the  self-consistent  plasma  electric  po¬ 
tential.  We  intend  to  demonstrate  that  this  new  simulation  tool  is  a  very  flexible 
and  efficient  alternative  to  the  Fokker-Planck  codes  that  are  very  widely  used  in 
the  plasma  physics  community.  Our  simulation  model  is  particularly  well  suited  to 
many  magnetic  fusion  configurations  (tandem  mirrors,  theta  pinches,  Q-machines, 
tokamak  divertors)  and  magnetospheric  double-layer  phenomena.  Our  first  appli¬ 
cations  will  be  to  the  study  of  electrostatic  particle  confinement  and  the  formation 
of  electrostatic  thermal  barriers  in  tandem  mirrors.  Figure  2  shows  the  particle 
and  energy  sources  and  sinks  in  one-half  of  a  symmetric  tandem  mirror  experiment. 
This  is  the  calculational  domain  of  our  axial,  one-spatial-dimension  model. 

Progress 

Significant  progress  has  been  made  in  development  and  testing  of  the  TESS 
(Tandem  Experiment  Simulation  Studies)  code  in  the  areas  of  code/ algorithm  op¬ 
timization  and  rudimentary  physics  applications.2  A  study  was  completed  in  which 
the  performance  of  two  alternative  implicit  differencing  schemes  was  compared  to 
determine  which  scheme  minimized  numerical  heating  and  cooling.  This  study  was 
based  upop  the  free  expansion  of  an  overall  charge-neutral  plasma  slab.  We  in¬ 
vestigated  operating  regimes  for  the  two  algorithms  with  respect  to  the  important 
parameters  upeAt  and  AxJXpe  »  where  upe  is  the  electron  plasma  frequency,  At  is 
the  time  step,  Ax  is  the  grid  spacing  and  Xpe  is  the  electron  Debye  length.  An  op¬ 
timum  value  of  Ax/\De  ss  3upeAt  for  upeAt  >  1  was  determined.  Values  of  u>peA t 
as  large  as  50  were  successfully  tested,  and  larger  values  can  be  accommodated. 
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Figure  2.  Schematic  drawing  of  the  1-D  axial  model 
ef  thermal  harrier  fermatioo  la  a  tandem 
mirror.  Indicating  the  principal  source* 
of  particles  and  heating. 


We  also  completed  and  successfully  tested  a  subroutine  for  binary  particle 
Coulomb  collisons.  This  major  physics  package  models  screened  Coulomb  collisions 
self-consistently.  i.  e.  the  distribution  function  evolves  due  to  the  collisions,  both 
momentum  and  energy  are  conserved,  and  the  entropy  of  the  plasma  either  in¬ 
creases  or  remains  constant.  Furthermore,  the  collision  kinematics  are  relativistic, 
and  the  Monte  Carlo  scattering  is  algorithm  is  vectorized.  A  complete  set  of  text¬ 
book  examples  of  collisional  relaxation  processes  have  been  simulated  and  reported. 
Detailed  descriptions  of  the  implicit  differencing  scheme  performance  study  and  the 
becnchmarking  of  the  binary  particle  collision  routine  will  be  presented  in  future 
Quarterly  Progress  Reports. 

Our  first  physics  applications  including  self-consistent  electric  potentials  have 
been  the  simulation  of  a  plasma  in  contact  with  tin  insulated  conducting  wall.  The 
observed  potential  drop  to  the  wall  and  its  dependence  on  electron  and  ion  tem¬ 
peratures  agree  with  analytical  sheath  models.  We  have  made  similar  simulations 
including  the  influence  of  applied  magnetic  mirror  fields  on  the  plasma  confine¬ 
ment.  These  studies  have  provided  a  very  valuable  testbed:  they  quantitatively 
establish  the  sensitivity  of  the  simulation  results  to  numerical  discretization  and 
statistical  noise,  and  clearly  demonstrate  that  TESS  successfully  combines  particle 
simulation,  Monte  Carlo  and  Fokker-Planck  simulation  tools. 

Future  Development 

The  implementation  and  testing  of  the  remaining  physics  packages  (wave  heat¬ 
ing,  neutral  beam  ionization  and  charge  exchanges,  axial  and  radial  particle  losses) 
over  the  next  several  months  will  enable  us  to  perform  realistic  thermal  barrier 


tandem  mirror  simulations,  thereby  establishing  credibility  in  the  physics  value  of 
TESS.  We  will  then  compare  TESS  to  existing  Monte  Carlo  and  Fokker-Planck 
codes.  WTe  will  also  continue  our  study  of  simulation  performance  characteristics 
and  explore  possibilities  for  multi-stepping  (sub-cycling)  and  orbit  averaging  of  the 
particle  data  to  improve  the  efficiency  of  the  code.3  Longer  term  goals  include  the 
extension  of  TESS  to  two  spatial  dimensions  and  different  magnetic  field  geome¬ 
tries. 

References 

1  A.  Friedman,  A.  B.  Langdon  and  B.  I.  Cohen,  ”A  Direct  Method  for  Implicit 
Particle- In- Cell  Simulation”,  Comments  Plasma  Phys.  Contr.  Fusion,  6, 
225(1981). 

2’  R.  J.  Procassini,  C.  K.  Birdsall  and  B.  I.  Cohen,  "Direct  Implicit  Particle 
Simulation  of  Tandem  Mirrors” .  Bull.  Amcr.  Phys.  Soc .,  31.  1621(1986). 

3  B.  I.  Cohen,  "  Orbit  Averaging  and  Subcycling  in  Particle  Simulation  of  Plas¬ 
mas”,  in  Multiple  Time  Scales (eds.  J.  U.  Brackbill  and  B.  I.  Cohen), 
Computational  Techniques,  Academic  Press,  London(1985). 


D .  Dynamic  Dimensioning  with  Heap  Management 
Lou  Ann  Schwager 

L  Introduction 

Dynamic  dimensioning  is  useful  for  those  who  need  to  automatically  increase  array  lengths  in 
a  computer  program  during  run  time.  In  particular  we  present  this  guide  to  demonstrate  the  use 
of  dynamic  dimensioning  in  particle  simulation  codes.  This  technique  allows  for  adding  memory  to 
arrays  just  before  their  pre-set  array  length  is  exceeded  during  a  run.  We  highly  recommend  that  the 
interested  user  acquire  a  copy  of  the  memory  management  section  in  the  BASCLIB  document  before 
attempting  any  of  this.  (Use  DOCUMENT  VIEW  BASELIB  with  keyword  “heap-management”  and 
keywords  for  individual  routines.) 

For  particle  simulation  codes,  certain  arrays,  such  as  position  and  velocity  for  each  particle, 
have  a  length  dependent  on  the  maximum  number  of  particles  npmax.  Usually  the  user  specifies 
npmax  in  the  code  prior  to  compilation  which  sets  the  lengths  of  the  arrays.  If  that  value  is  too 
small,  then  the  run  terminates  (or  generates  an  error  downstream)  when  array  lengths  are  exceeded. 
Then  the  entire  computation  costs  are  wasted.  Also  for  a  simulation  model  that  generates  new 
particles  with  time  and  starts  void  of  particles  then,  through  most  of  the  run,  much  unused  space 
exists  in  the  particle  dependent  arrays.  Those  users  familiar  with  the  NMFECC  system  realise 
that  a  larger  program  field  length  requires  either  a  higher  priority  (and  thus  a  higher  cost  charged 
for  computation)  or  more  personal  time  invested.  To  mitigate  these  two  problems,  using  dynamic 
dimensioning  can  avoid  the  case  of  premature  termination  and  can  reduce  the  average  cost  of  a 
run.  Especially  improved  are  those  runs  requiring  over  10,000  total  particles  and  over  ten  minutes 
of  Cray-1  computation  time.  Memory  management  does  cost  some  i/o  and  system  time  so  that 
optimally  it  should  be  called  on  only  a  few  times  during  a  run. 

II.  Technique  Implementation 

Usually  array  lengths  are  specified  in  the  code  prior  to  compilation.  Hence  to  use  memory 
management,  before  compilation  we  set  the  length  of  these  expandable  arrays  to  one.  In  the  coding 
before  any  of  these  arrays  are  referenced,  the  initial  values  for  these  lengths  are  read  in  from  the 
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input  file.  Then  we  call  on  the  heap  manager  to  allocate  a  moveable  block  of  memory  for  each 
expandable  array. 


Before  embarking  on  how  to  accomplish  the  initial  allocation,  we  must  explain  the  use  of  pointers 
in  memory  management.  (Pointers  are  discussed  in  the  CFT  document  under  keyword  ‘‘pointer- 
statement’’.)  Pointers,  which  have  integer  type,  provide  the  base  address  of  an  array  or  variable. 
Changing  the  value  of  a  pointer  moves  the  location  of  the  entire  array.  The  pointer  statement  for 
a  particular  array  must  reside  in  the  routine  of  use.  This  restriction  exists  because  the  pointered 
array  cannot  be  passed  into  a  routine  as  a  dummy  argument  or  in  a  common  block;  whereas  the 
pointer  itself  can  be. 

Pointers  can  be  passed  through  labeled  common  blocks  to  the  subroutine  where  memory  is 
managed.  Using  the  BASELIB  routine  MZMALLOC  for  each  array  assigns  the  initial  moveable 
block  of  memory.  The  first  argument  of  MZMALLOC  is  the  pointer  of  the  array  and  the  second 
argument  is  the  block  length  in  words.  The  array  can  alao  be  dimensioned  in  the  pointer  statement. 
Note  that  error  status  is  returned  when  almost  any  of  the  MZ  routines  are  called.  When  the  error 
status  equals  sero,  the  task  has  completed  successfully.  Bence  nonsero  values  should  be  trapped 
and  the  run  suspended  or  terminated.  (See  the  BASELIB  document  for  details.)  At  this  point  in 
implementing  the  above  changes,  the  user  can  test  the  program  by  choosing  a  value  of  npmax,  in 
the  input  file,  which  will  not  be  exceeded  during  a  run. 

As  the  run  proceeds  and  array  values  are  assigned,  a  test  a— cases  whether  the  array  length  is 
about  to  be  exceeded.  When  this  occurs,  that  routine  should  call  the  subroutine  which  lengthens 
the  array  blocks.  (Placing  this  expansion  routine  in  a  separate  subroutine,  though  not  required, 
facilitates  debugging.)  The  BASELIB  routine  MZCHANGE  expands  an  allocated  block  of  memory. 
Its  first  argument  is  the  pointer  for  the  array  and  the  second  is  the  total  new  block  length.  The 
BASELIB  routine  MZSTATUS  can  be  called  on  to  provide,  for  curiosity’s  sake,  the  free  space 
available  after  each  MZCHANGE  call.  The  BASELIB  routine  IZPROGFL  provides  the  program 
field  length. 
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In  printing  all  the  pointer  values,  program  field  length,  and  available  free  space,  the  curious 
user  can  observe  how  the  heap  manager  shuffles  the  newly  expanded  arrays.  Usually  with  the  first 
call  to  the  expansion  subroutine,  the  first  array  changed  is  moved  to  an  address  larger  than  that  of 
the  last  array.  (“Last”  is  defined  as  the  final  array  initially  allocated  memory  with  MZM ALLOC.) 
The  second  array  changed  may  move  beyond  the  first  array  or  may  move  into  the  old  space  of 
the  first  array.  This  shuffling  depends  on  the  amount  of  extra  length  added  and  the  array  site. 
Obviously  after  several  calls  to  the  expansion  subroutine,  the  relative  base  addresses  of  these  arrays 
are  considerably  scrambled  from  the  original  positions.  Consequently  if  the  user  employs  a  pointer 
which  depends  on  the  value  of  an  array’s  base  address,  then  this  pointer  must  be  updated  after  all 
changes.  The  heap  manager  only  returns  the  new  value  of  the  array’s  base  address. 

We  observe  that,  with  only  using  MZCHANGE,  after  the  heap  manager  juggles  all  the  arrays 
in  the  expansion  routine,  a  considerable  amount  of  free  (unallocated)  space  exists.  The  BASELIB 
routines  MZPARAM  and  MZPACKH  are  used  to  reduce  this  free  space.  Thus  the  minimum  program 
field  length  is  maintained  during  the  run  until  the  expansion  routine  may  again  be  called  Deciding 
how  and  when  to  use  these  two  routines  requires  some  operating  experience. 

The  user  sets  the  second  argument,  epsilon,  and  the  fourth  argument,  maxfat,  in  the  array 
required  as  input  for  an  MZPARAM  call.  This  epsilon  value  sets  the  amount  of  extra  space  allowed 
to  accompany  a  changing  block.  Consultants  recommend  that  it  be  set  to  no  less  than  the  shortest 
block  length.  The  value  maxfat  specifies  the  maximum  space  allowed  to  accumulate  at  the  end 
of  the  heap  before  the  heap  manager  removes  the  excess  and  shortens  the  program  field  length. 
The  default  for  maxfat  is  32.  Because  the  heap  manager  can  shorten  the  field  length  without 
provocation,  we  fix  maxfat  initially  to  a  large  value  such  as  the  largest  block  sise.  Then  after  all 
arrays  are  expanded  (at  the  end  of  the  expansion  routine)  we  set  maxfat  back  to  32  and  then  call 
MZPACKH 

Next  calling  MZPACKH  repacks  the  heap  to  move  free  space,  according  to  epsilon,  to  the  end 
of  the  heap.  Then  the  heap  manager  automatically  removes  the  extra  end  space  (beyond  maxfat). 
The  argument  for  MZPACKH  is  set  equal  to  sero.  As  it  currently  stands,  MZPACKH  will  shuffle  a 
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block  of  memory  (to  consolidate  free  apace)  only  if  sufficient,  contiguous  free  space  exists  to  shuffle 
the  allocated  block  into.  In  other  words,  MZPACKH  is  not  smart  enough  (as  yet)  to  slide  blocks 
around  the  small  extra  spaces,  even  if  epsilon  were  set  less  than  the  smallest  block  length 

With  dynamic  dimensioning,  we  provide  the  caveat  that,  although  array  space  may  be  allo¬ 
cated,  strange  errors  ensue  if  an  array  element  is  referenced  before  its  value  is  initialised  Because 
the  garbage  residing  at  the  referenced,  uninitalised  location  may  propagate  through  the  code,  an 
overflow/underflow  error  may  appear  several  subroutines  later. 

III.  Technique  Demonstration 

The  particle  simulation  code  PDW1  is  used  to  illustrate  the  dynamic  dimensioning  technique. 
In  PDW1  the  position  and  velocity  arrays  are  sssigned  pointers  in  the  routine  MAIN .  Other  particle- 
dependent  arrays,  such  as  flags  for  certain  locations  within  the  simulation  system,  are  pointered  in 
the  subroutine  where  these  flags  an  set.  For  example,  in  MAIN  the  pointer  statement  for  paeition 
array  x  is  pointer(px,  x(  1 )) .  The  pointer  px  passes  through  a  common  block  to  subroutine  START 
where  the  heap  manager  allocates  memory  via  ierr  as  mzmalloc{px ,  npmax  •  nsp)  where  nsp  is  the 
number  of  species.  Memory  for  ail  the  expandable  arrays  is  allocated  initially  in  subroutine  START 
In  PDW1  the  test  for  exceeding  the  maximum  number  of  particles  npmax  occurs  while  particles 
are  added  during  a  call  to  LOAD  or  ADJUST.  If  the  test  is  true,  then  subroutine  EXPAND  is  called. 

Suppose  that  the  particle-dependent  arrays  an  position  x,  velocity  v,  and  some  flag  if  The 
x  and  v  arrays  have  length  npmax  *  nsp  and  if  has  length  npmax.  The  x  and  v  arrays  contain 
the  values  for  species  stacked  in  one-dimenaioa  (rather  than  in  nep-dimensions).  Consequently  the 
elements  in  array  nloc.  which  mark  the  beginning  of  each  species,  must  also  be  increased  with  the 
new  npmax  after  EXPAND  is  called.  The  expansion  routine  EXPAND  follows. 
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subroutine  expand(npmax,nsp,nloc) 
use  maxpar 

c*  Below  are  the  three  pointers  which  locate  the  base  address 
c*  of  each  respective  array:  *,  v,  and  if 
common  /pmain/  px,  pv 
common  /pfla g/  pif 
common  /shuffle/  nadd 
dimension  nloc(nsmaxl),  kparam(4) 
data  kparam(l)/-l/,  kparam(3)/-l/,  kpackh/0/ 
integer  px,  pv,  pif 

c*  JV add,  the  memory  added,  is  the  user’s  best  guess  which 
c*  affects  the  number  of  times  EXPAND  may  be  called 
nadd  =  5000 

c*  Update  npmax  and  nloe 
npmax  =  npmax  +  nadd 
do  10  i=l,nsp 
nloc(i+l)  -  1  *  npmax 
10  continue 

c*  Set  max  fat  and  eptilon 
kparam(2)  =  npmax 
kparam(4)  =  npmax  *  nsp 
iparam  *  msparam(kparam,4) 
c*  Expand  each  array  and  print  status, 
c*  The  error  value  ichanft  is  incremented  by  ten  with 
c*  each  rnxchange  to  aid  in  locating  errors  and  in 
c*  understanding  the  status  printout. 

e*  With  a  correctly  completed  task,  the  mzchangc  value  is  sero. 
•change  =  mschange(px,npmax*nap)  +  10 
if(ichange.gt  10)  goto  20 
call  statusfichange) 

■change  =  mschange(pv,npmax*nsp)  +  20 
iffichange.gt  20)  goto  20 
call  status!  ich an ge) 

•change  =  rnschange<pif,npmax*nsp)  +  30 


if(ichange.gt.30)  goto  20 
c*  Reset  maxfat  to  default  and  repack  heap 
kparam(4)  =  32 
iparam  =  miparam(kparam,4) 
ipackh  =  mzpackh(kpackh) 
call  status(ichange) 
go  to  30 
20  call  finish 
30  continue 
return 
end 

subroutine  status(ichange) 

c*  Print  program  length,  free  space,  and  pointer  values 
c*  in  OUTPUT. 

common  /pmain/  px,  pv 
common  /pflag/  pif 
dimension  aprog(3),  kstat(2) 
iprog  =  izprogfl(aprog) 
is  tat  =  msstat(kstat,2) 

write(3,ll)  kstet(2),eprog(2),px,pv,pif,iebange 
11  format( lx, “STATUS:  free  spaces”  ,i6, “field  lengths”  ,i8,/, 

*  “pxs”  ,i8,“pv=”  ,i8,“pif=”  ,i8,“for  ichanges”  ,i3) 
return 
end 

With  x  and  v  newly  expanded,  if  the  program  adds  a  particle  of  the  first  species,  then  the  first 
element  in  x  and  v  for  the  second  species  will  be  overwritten.  Consequently  we  must  distribute  the 
extra  space  allocated,  which  is  now  residing  at  the  end  of  each  array,  internally  to  the  end  of  each 
species.  Immediately  after  calling  EXPAND,  we  then  call  the  subroutine  SHUFFLE  which  follows 
below. 
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subroutine  shuf9e(s,nloc,nsp) 
c*  Redistribute  memory  for  species-dependent  array  a 
common  /shuffle/  nadd 
dimension  s(l),  nloc(l) 

c*  Routine  MOVE  (in  FORTLIB)  moves  array  sections  to 
c*  new  addresses  within  that  array. 

c*  Note  that  in  PDW1  the  suboutine  MOVE  must  be  renamed, 
imove  =  nloc(2)  -  nadd 
do  10  isp=nsp,2,-l 
moveto  =  nloc(isp)  +  1 
movefrm  =  nloc(isp)  -  (nadd*(isp-l))  +  1 
call  move(s( moveto) ,s( movefrm), imove) 

10  continue 
return 
end 

For  any  species-dependent  array,  subroutine  SHUFFLE  must  be  called  just  prior  to  its  next 
reference  during  the  run.  For  example,  if  such  an  array  is  in  subroutine  OUTPUT,  then  it  is  not 
referenced  right  after  the  call  to  EXPAND  (in  ADJUST  or  LOAD)  but  perhaps  many  time  steps 
later.  Hence  we  can  set  a  flag  in  EXPAND  which  is  passed  via  a  named  common  block  to  OUTPUT. 
When  this  flag  is  true  then  we  call  SHUFFLE  (from  OUTPUT)  and  then  reset  the  flag  to  false. 
IV.  Conclusions 

Hopefully  the  above  discussion  covers  most  of  what  the  user  may  need  to  undertake  the  use 
of  dynamic  dimensioning.  This  technique  can  also  be  used,  for  example  in  matrix  inversion,  to 
allocate  and  deallocate  memory  should  a  large  array  be  temporarily  required.  The  whole  process 
also  educates  the  user,  who  previously  was  familiar  only  with  the  old  linkage  (link=0)  system,  on 
the  strange  new  creature  of  the  link=2  system  called  the  heap  and  its  management. 
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